The impact of dietary calcium and phosphorus on mitochondrial-linked gene expression in five tissues of laying hens

Mitochondria and the energy metabolism are linked to both, the availability of Ca and P to provide the eukaryotic cell with energy. Both minerals are commonly used supplements in the feed of laying hens but little is known about the relationship between the feed content, energy metabolism and genetic background. In this study, we provide a large-scaled gene expression analysis of 31 mitochondrial and nuclear encoded genes in 80 laying hens in the context of dietary P and Ca concentrations. The setup included five tissues and gene expression was analysed under four different diets of recommended and reduced Ca and P concentrations. Our study shows, that mitochondrial gene expression is reacting to a reduction in P and that an imbalance of the nutrients has a higher impact than a combined reduction. The results suggest, that both strains (Lohmann Brown and Lohmann Selected Leghorn) react in a similar way to the changes and that a reduction of both nutrients might be possible without crucial influence on the animals’ health or gene expression.


Introduction
Phosphorus (P) is an essential mineral for all living organisms which must be continuously supplied and is needed in poultry for growth, health and the energy metabolism [1]. The ability of poultry to degrade the natural present phytic-acid (InsP 6 ) is limited [2,3] and thus feed supplements derived from rock phosphate are added to maintain the P supply. The availability of rock phosphate is limited [4] and thus a reduction of its usage is of utmost interest.
Another important mineral essential for laying hens is calcium (Ca) [5] which is needed e.g. to form eggshells. Recent studies in broiler chickens have shown, that endogenous InsP 6 degradation is reduced when mineral P and Ca are supplemented [6][7][8] and there is a strong interaction of P and Ca content regarding egg-shell quality and the number of produced eggs [9]. Many studies suggested, that the recommended dietary P content in the feed of laying hens might be too high, and can be reduced, without significant negative effects on performance and health of the animals [5,[10][11][12]. Hence, a better understanding of the effects of dietary P and Ca is necessary to implement these suggestions and adjust their dietary provision accordingly.

PLOS ONE
PLOS ONE | https://doi.org/10.1371/journal.pone.0270550 June 24, 2022 1 / 17 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Both used strains (Lohmann Brown and Lohmann Selected Leghorn) are commercial important and selected for egg production [13][14][15]. Despite their similarities in egg production previous studies have shown differences between them concerning body weight, gene expression and phytate degradation in the context of the productive life span and changes in dietary Ca and P concentrations [10,16,17]. These results suggest differences in the strains reaction to dietary changes, which will be analysed in this work.
In this study we focus on mitochondrial gene expression since mitochondria are linked to P and Ca availability as well as to the animals' fitness and energy metabolism. Mitochondria are the main energy producers of the cell through the process of oxidative phosphorylation (OXPHOS) [18] and this process depends on the availability [1] and is influenced by the concentration [19] of P. In addition, they play a major role in regulating Ca 2+ homeostasis [20], which is an important factor in cell signalling since Ca 2+ controls many cellular functions [21] including gene expression [22] and the regulation of OXPHOS [23].
In our experimental setup, we included all mitochondrial encoded OXPHOS subunits, as well as nuclear encoded ones such as NDUFB6, which is a subunit essential for the electron transport in complex I of the respiration chain [24] and SOD2, which detoxifies reactive oxygen species (ROS) produced during OXPHOS in mitochondria [25,26]. In addition, we included nuclear genes which are part of the regulatory network, linking mitochondrial gene expression and biogenesis to external stimuli, such as PGC1α [27], nutrient sensitive factors such as MTOR as well as subunits of AMPK, which is part of the adaptive response to energy deficit [27]. Another important key player is IGF-1α, which has been described as participating in P transport [28], and a reduced IGF1 expression increase the effect of Ca deficiency on bone accretion in mice [29]. IGF-1α has also been linked to body weight in chickens [30], which makes it a promising gene in our experimental setup.
Preliminary studies mostly focus on physiological traits important for agricultural purposes, but also try to understand the mechanisms behind the effects of dietary P and Ca. In a study focusing on phytate degradation, transcellular mineral transporters, and mineral utilization of the same hens [10] no P-mediated effects were identified and a major question rising from this studies is, how the animals react to compensate the reduced amount of P and Ca. In this study, we focus on the expression of mitochondrial-linked and nutrient sensitive genes and test the following hypotheses: 1. Mitochondrial gene expression reacts to the reduction of P and Ca 2. The strains react differently to the changes in the diet composition 3. Genes regulating mitochondrial gene expression and biogenesis react to the dietary changes Our study benefits from the already published analysis of phenotypic traits [10] as well from the known genetic background of the two strains [31]. Together with the possibility of a controlled environment during the experimental procedure it was possible to detect and analyse gene expression changes in the context of dietary adjustment.

Animals and experimental setup
The animal experiments were performed at the Agricultural Experiment Station of the University of Hohenheim, Germany. They were approved by the Regierungspräsidium Tübingen, Germany (Project no. HOH50/17TE) in accordance with the German Animal Welfare Legislation.
We used 80 laying hens: 40 Lohmann brown classic (LB) and 40 Lohmann LSL-classic (LSL) white leghorn hybrids contributed by Lohmann Tierzucht (Cuxhaven, Germany). The hens originated from an experiment addressing the utilization of phosphorus (P) and calcium (Ca) under different dietary conditions [10]. The experimental setup is described in detail in Sommerfeld, Omotoso, et al., 2020 [10] and will only be outlined briefly in the following.
The hens were reared together under standard conditions, with diets according to the requirements of each period, based on soy and corn meal. Ten father lines per strain were selected prior to the start of the experimental phase. After 27 weeks four hens per rooster were chosen and placed individually in metabolism units (1m × 1m × 1m) where the hens received specific diets for the following three weeks. Four different feed compositions were used (Table 1) and fed ad libitum, each group contained one hen per father line.

Samples and RNA extraction
Samples of five tissues (breast muscle, ileum, duodenum, liver and ovary) were taken on four consecutive days with random distribution of the four diets in week 31. The animals were individually stunned with a gas mixture of 35% CO 2 , 35% N 2 , and 30% O 2 and killed by decapitation at the Agricultural Experiment Station of the University of Hohenheim [10]. The samples were directly taken after slaughtering and were immediately placed on dry ice and stored at -80˚C until RNA extraction.
RNA was extracted using TRIzol Reagent (Thermo Fisher scientific Inc., Massachusetts, USA) according to the manufacturers' instructions with modifications described in Dreyling and Hasselmann 2022 [32]. Samples were dissolved in nuclease-free water, RNA concentration and quality in form of 260/280 and 260/230 ratios were measured using a NanoDrop 2000/ 2000c Spectrophotometer (Thermo Fisher scientific Inc., Massachusetts, USA). In addition to the provided NanoDrop values for all samples, a representative subset including samples of different quality and quantity was measured on a Qubit 4 (Thermo Fisher scientific Inc., Massachusetts, USA) using the Qubit RNA IQ Assay Kit (Thermo Fisher scientific Inc., Massachusetts, USA). The samples were stored at -80˚C until further processing.

Real time PCR
The primer design and assay evaluation are already published in Dreyling and Hasselmann (2022) [32], and were not repeated specific for this experiment. The experimental procedure is identical and only described in brief in the following. In this study, 28 candidate genes were used, whereas now five additional genes have been integrated into the final set of primers: ATPF0, ND2, ND3, PRKAB1, and PRKAG3. A list with the gene names and their abbreviations used in this study can be found in Table 2.The used primer set including product size, primer efficiency and accession numbers of the reference sequences can be found in S1 Table. In addition, melting and standard curves are provided which were produced the same way as described in Dreyling and Hasselmann (2022) [32]. All analyses were performed on a Biomark HD system (Fluidigm Corporation, San Francisco, USA), following the protocols for gene expression analysis using five 96.96 IFCs using the Delta Gene Assays protocol with the manufacturers standard protocol for fast PCR and melting curve as described in S2 Table. The protocol includes a DNase digestion prior to the reverse transcription and a pre-amplification with multiplexed primers followed by an Exonuclease I treatment. Pre-amplification bias was already described and discussed in Dreyling and Hasselmann (2022) [32]. All qPCRs were performed in duplicates and the samples were placed randomly on the chips, only grouped by individual to avoid any bias of sample arrangement. An internal control and a negative control (throughout all preparation steps) were included to detect variance between the runs and potential contamination. In each well of the IFC 2.25μl of the diluted Exo I digested sample were added, resulting in 3.015nl in each reaction chamber. Information about the primer-pairs and thermal cycling conditions can be found in Table A-C in S1 File.

Data preparation
Quality control. For data evaluation and quality control the Fluidigm Real-Time PCR analysis software was used. Only Cq-values from reactions with logarithmic increase of fluorescence and specific melting points were used for the following analyses. After the automatic quality check of the software, the results were evaluated by eye and revised manually if necessary. The quality threshold was set to 0.65 and the peak ratio threshold to 0.8.
The results of the internal control were checked to detect possible variation due to technical issues.
Reference gene evaluation. We used PPIA and ACTB as reference genes for normalization. Our previous evaluation already showed that GAPDH is strongly influenced by tissue and thus we included it as candidate gene in our study. To verify that GAPDH shows high variance between tissues in this experiment as well, we performed a reference gene evaluation using Normfinder [33]. As input one individual of each strain and diet was used, including all five tissues per individual to cover all our variables of interest. Normalization was tested for tissue type and diet.
Calculating relative gene expression. Means of duplicates were calculated of all samples with two successful runs. For samples that only had one successful duplicate this run was used.
Δct was calculated as the difference between the average cycle threshold (ct) of the internal control to the ct of the corresponding sample. E refers to the converted primer efficiency and GOI to gene of interest. RQ are relative quantity values calculated using E and Δct. PPIA and ACTB were used as reference genes (refs). The data was log2 transformed prior to statistical analysis.

Statistical analysis
Linear mixed models were implemented and adjusted to each of the genes.
Where Y is relative gene expression, ε is the residual error; strain, diet and tissue are fixed effects, with individual and father as random effects. All modelling was performed in R (R Core Team 2019, Version 3.6.1) using the lmerTest [37]. A three factorial analysis of variance (anova) was performed to evaluate the influence of fixed effects and pairwise tukey posthoc tests (package emmeans [38]) were performed to detect differences between strain, diet and tissues in various combinations. The output of the model is estimated marginal means (emmeans), which are used for statistical analyses throughout the whole study. The fulfilling of normal distribution and the homogeneity of variance were evaluated using QQ and residual plots. Outliers were removed for each dataset using the interquartile range, except for SDHA, PRKAA2, PRKAB2 and GAPDH because for these datasets a removal of outliers would have included too many samples to perform a proper analysis.

Results
After the quality filtering, we received a dataset of 13,487 ct values from 35 Genes and all 400 samples. From 184 samples we obtained high-quality ct values for all genes, whereas for the remaining samples at least one gene was missing. In the final analyses all samples were included since no relationship between sample quality and the successful run of all genes was recognizable. Tables containing concentrations and quality of the RNA extracts can be found in S3 File.
The statistical model revealed, that no included gene was significantly influenced by the strain, two genes were influenced by diet and all genes by the tissue. The most frequent interaction was between strain:tissue (6), followed by strain:diet and strain:tissue:diet (1) while there was no significant interaction of tissue:diet. The results of the tests derived from the statistical models for all genes and factors can be found in S1 Table. Sample numbers per gene, strain, tissue and diet can be found in S2 and S3 Tables.

Strain differences
The two-way hierarchical clustering analyses revealed a different number of clusters within each strain of laying hens. Two clusters were calculated by the cubic clustering criterion [39] for both strains (Fig 1). The analysis was performed on 92 samples of each strain, including 19 to 25 samples per diet per strain and 14 to 24 samples per tissue and strain. In general, both strains showed differentiation different tissue types, especially breast and liver, with more pronounced tissue specific expression pattern in the LSL strain. Consistent for both strains, a strikingly high up-regulation of three subunits of AMPK (PRKAA2, PRKAB2 and PRKAG3) and GAPDH are found within breast tissue. The number of clusters was two in the LB and four in the LSL strain, which reflects the clearer separation of tissue types in the LSL strain.
Using our statistical linear mixed model, we tested for the overall impacts of strain, tissue and diet on gene expression. We observed no difference between the two strains for all genes under all four nutritional conditions, except for IGF-1α and UQCRC1 (Fig 2). Both genes were showing higher gene expression in the LSL strain for all four diets, with a significant difference under low P and high Ca (diet 4) for IGF-1α (p = 0.0076), and low Ca and low P conditions (diet 2) for UQCRC1 (p = 0.0214).
For IGF-1α these strain differences are most pronounced in three tissues: breast (higher in LB hens under low P conditions (p = 0.0227 for P-Ca-and p = 0.0415 for P-Ca+)), liver (always higher in the LSL strain, p = 0.003 for P+Ca+, p = 0.0001 for P-Ca-and p<0.0001 for P+Caand P-Ca+) and ovary tissue with contrasting pattern among the diets under P+Ca-(p = 0.0116) and P-Ca+ (p = 0.0037) (Fig 3). In ovary tissue, the significant decrease of gene expression in the LB strain under P-Ca+ is consistent with the overall trend (Fig 5) while the increase in the LSL strain does not follow this pattern.

Differences between the diets
The overall gene expression (mean of all genes, tissues and strains per diet) was lowest under low P and high Ca conditions; however, the differences between the diets were not significant. The statistical model revealed, that the expression of only two genes (SOD2 and NDUFB6) was affected by the diet. Nevertheless, with pairwise comparisons of the diets, more genes showed significant gene expression differences: ND3, CytB, SOD2, COXC6, and NDUFB6. For all these genes, the expression was lowest under low P and high Ca conditions. The expression of SOD2, COXC6 and NDUFB6 was significantly higher under high P compared to low P under high Ca conditions (Fig 4) for both strains analysed together. In addition, the LB strain showed expression differences in CytB, where the expression was higher under high P and Ca levels as well as under low P and Ca levels compared to diet 4 (P-Ca+).
In addition to the already described differences between the strains in the expression of IGF-1α, the differences between P+Ca-and P-Ca+ in ovary tissue were significant in the LB strain (p = 0.0062) and strong but not significant in the LSL strain (p = 0.0563) (Fig 3). Additionally, the expression in the LB strain was significantly higher in liver tissue under P-Cacompared to P-Ca+ (p = 0.0088)

Tissue differences
As already discussed in the context of time dependent gene expression [32], the tissue had the strongest influence in our experimental setup, influencing the expression of all included genes. The gene expression was significantly highest for all genes in breast muscle compared to the remaining four tissues (p<0.001 for all except PRKAA1 where p = 0.048 when breast compared to liver), except for IGF-1α and PRKAG2, where there was no difference compared to  PLOS ONE liver tissue. The expression was lowest in all tissues when fed diet 4 (P-Ca+) (Fig 5). The already described diet-dependent gene expression changes in SOD2 and NDUFB6 were only observed in breast muscle tissue, where the expression was significantly lower under P-Ca + than under P+Ca+ levels (for more details see S3 Table). In GAPDH highly significant differences between all pairwise tissue comparisons were observed (p<0.0001), which supports the decision to abandon it as a reference gene in our study.

Discussion
In our experimental setup we were able to analyse a vast number of mitochondrial and nutritional linked genes in the context of dietary changes in P and Ca contents. We hypothesized, that one compensatory mechanism of changes in the diets is the adaption of mitochondrial gene expression, since it is directly linked to the availability of P and the fitness of the individuals.

Differences between the strains
For none of our candidate genes, the gene expression differed between the strains, and only two genes (IGF-1α and UQCRC1) showed different gene expression in specific diets (under low P concentration). These observations indicate, that both strains react to the changes in dietary P and Ca content the same way and also in the hierarchical clustering analyses, the pattern of both strains was similar. A genome wide gene expression analysis comparing the same strains as included in this work identified genes related to the GO-cluster of phosphorous metabolism (GO-IDs: GO:006468, GO:0006793, GO:0006796, GO:0016310) to be down regulated in hens of the LSL strain [14] and Sommerfeld et al., 2020 [10] identified two sodium/

PLOS ONE
phosphate co transporters to be higher expressed in the LSL strain using the same individuals as in this work. These results indicate, that there are differences between the strains related to P metabolism, however the genes included in this work showed no general differences between the two strains. It must also be noted, that Sommerfeld et al. 2020 [10] identified no interaction of strain and diet on the included genes, which supports the hypothesis that the reduction of P and Ca content in the used diets was not sufficient to improve the expression levels of our gene targets significantly. Even if we were able to show differences in gene expression between strains and diets in some genes in our study, the majority of the included genes showed no reaction to the dietary changes. A significant change in the mineral concentrations with an detrimental effect on the whole animal might lead to stronger effects on the animals as described in this work or by Sommefeld et al. 2020 [10].

Mitochondrial gene expression in the context of dietary changes
Since both nutrients are linked to the mitochondrial energy-metabolism [1,23] we suggested an adaption of mitochondrial gene expression as a compensatory reaction to the changing amount of the minerals. Our data revealed, that most genes showed no significant difference in gene expression according to the changes in the diet. However, the gene expression was lowest under P-Ca+ conditions, while the reduction of both minerals or only Ca had a smaller effect on gene expression. These results suggest, that the effect of reduced P concentrations is stronger, when there is an imbalance of the proportion, especially under low P.
Four of the five genes showing significant differences were part of the electron respiratory chain, representing complexes I, III, and IV and the ROS detoxifying gene SOD2. The gene expression was significantly lower under P-Ca+ compared to P+Ca+ in genes representing subunits from Complex I and IV of the respiration chain and the ROS detoxifying gene SOD2. The same observation was made for CytB in the LB strain (Fig 4). This observation suggests that the reduced availability of P impacts distinct parts of the respiratory system, resulting in reduced gene expression. Since it is known, that the expression of whole complexes can be regulated by the expression of individual subunits [41,42], the reduced expression of single subunits might regulate the amount of the whole complex. Additionally, previous studies showed pattern of co-expression of the OXPHOS complexes [43], which is also shown in our analysis under low P and high Ca conditions. A general reduction of assembled OXPHOS complexes might be the result of the observed expression pattern. The potential reduced production of ROS resulting from diminished OXPHOS activity leads to a reduced need of SOD2. SOD2 expression has been linked to impact immunity against bacterial infections in zebrafish [44] and ROS are known to play a role in the reaction to inflammatory disease [45][46][47]. Thus, the differences in SOD2 expression in laying hens might indicate differences in resistance to infections as well. In addition, an increase in the production of ROS in the mitochondrion is linked to the process of ageing in many species [48] and the increase of SOD2-expression protects the mitochondrion from damage, which would otherwise lead to the death of the cell (as stated in Santos et al., 2018 [48], Yin et al., 2018 [49] and cited references within). Regarding the missing of differences in gene expression between the other treatments suggests that a reduction of P alone is more crucial than a reduction of P and Ca or Ca.

Gene expression differs between different tissue types
As described in the context of life span [32] the gene expression in breast muscle was significantly higher than in all other tissues, followed by liver tissue. Studies in humans and chimpanzees have shown, that the differences in gene expression are higher between tissues than between species and suggests to include different tissue types in functional genomics studies [50]. In liver, the high rates of gene expression might be explainable by the high amount of functions of the tissue, ranging from the conversion of glucose into glycogen, the filtration of blood to the production of cholesterol [51]. In the duodenum and ileum the mucosa and microbiome are important participants in the process of digestion and uptake of nutrients [52], which might explain the lower rates in the tissue itself. However, the high activity in breast muscle tissue is surprising, especially since the process of growing was already finished during sampling and the purpose of the strains focuses on egg laying instead of meat production. In general, there are not many studies comparing gene expression rates between tissues since most studies focus on differences between different treatments, diseases or developmental states.

IGF-1α in the context of nutritional changes
We observed most differences in gene expression between strains and diets in the growth factor IGF-1α. This gene plays a major role in a variety of tissues and functions, e.g. metabolic homeostasis [53,54], and growth [55,56] and nutrition has been identified as a key factor of IGF1 regulation in humans [57,58]. This sensitivity might be the reason of the significant differences in expression between the diets observed in this study, especially since malnutrition is known to reduce circulating IGF1 in mice [59]. Additionally, malnutrition has been linked to decreasing IGF1 expression in liver tissue [60], which is reflected by the significantly lower expression under P-Ca+ compared to P-Ca-in liver tissue of the LB strain in this study. This observation also leads to the conclusion, that an imbalance of the minerals is detrimental compared to the reduction of both minerals, and is thus a form of malnutrition. An interaction of P and Ca content have also been shown in the context od egg-shell quality and quantity of eggs [9]. Simultaneously the gene expression is significantly lower under P-Ca+ conditions in the LB strain, which is another indication of differences in the reaction to the reduction of Ca in both strains. The most prominent difference between the strains was the contrary expression of the strains under P+Ca-and P-Ca+ conditions. Even if it is long known, that IGF-1α plays a role in avian ovaries [61], we could not observe any changes in egg weight between the strains matching the change in expression pattern [10].

Expression changes in nutrient sensitive and mitochondrial regulatory genes
We included nutrient sensitive genes such as MTOR and AMPK and the nutrient sensitive mitochondrial regulator PGC1α in our study and the missing reaction to our dietary changes is striking. All of these genes have been analysed in the context of changes during the productive life span of laying hens [32] using the same technical approach, which makes it rather unlikely, that our setup is failing in detecting expression changes. The same hens have been analysed in the context of performance traits such as body weight, feed intake, average egg weight and P/Ca efficiency, where no diet specific changes could be observed [10]. In accordance to other studies [5,11,12] the authors conclude, that a 20% reduction of P is not affecting the animals and thus, the recommended concentrations in the feed of these animals might be too high. The missing effect on genes that are part of the mitochondrial regulatory network supports this hypothesis.

Conclusion
We performed a large-scaled analysis of mitochondria-linked gene expression in laying hens in the context of P and Ca content in the diets. Our study revealed interesting differences of gene expression of subunits covering most OXPHOS complexes under low P and normal Ca concentrations in the diets. Together with the decrease in the expression of the ROS detoxifying gene SOD2, an interesting part of the regulation of mitochondrial gene expression has been revealed. In addition, the effects on the growth factor IGF-1α showed that the reduction of P in the diet has an effect on the mitochondrial regulatory network as well. We also observed that an imbalance of both minerals seems to have greater influence of gene expression than the reduction of both nutrients, especially under low P conditions.